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Abstract 

We introduce a new class of arbitrary-order exponential time differencing methods based on 
spectral deferred correction (ETDSDC) and describe a simple procedure for initializing the requisite 
matrix functions. We compare the stability and accuracy properties of our ETDSDC methods to 
those of an existing implicit-explicit spectral deferred correction scheme (IMEXSDC). We find that 
ETDSDC methods have larger accuracy regions and comparable stability regions. We conduct 
numerical experiments to compare ETD and IMEX spectral deferred correction schemes against a 
competing fourth-order ETD Runge-Kutta scheme. We find that high-order ETDSDC schemes are 
the most efficient in terms of function evaluations and overall speed when solving partial differential 
equations to high accuracy. Our results suggest that high-order ETDSDC schemes are well-suited 
to work in conjunction with spectral spatial methods or other high-order spatial discritizations. 
Additionally, ETDSDC schemes appear to be immune to severe order reduction, a problem which 
affects other ETD and IMEX schemes, including IMEXSDC. 

Keywords: Spectral deferred correction, exponential time differencing, implicit-explicit, high-order, 
stiff-systems, spectral methods. 


1 Introduction 

In this paper we present a new class of arbitrary-order exponential time differencing (ETD) methods 
for solving nonlinear evolution equations of the form 

= K(j) + A/'(t, (j)) 

where A is a stiff linear operator and A/" is a nonlinear operator. Such systems commonly arise when 
discretizing nonlinear wave equations including Burgers’, nonlinear Schrodinger, Korteweg-de Vries, 
Kuramoto, Navier-Stokes, and the quasigeostrophic equation. ETD Adams methods O [20], ETD 
Runge-Kutta methods (T] [22l [25l [181 El ESI ED] , and ETD general linear methods |34l [20] are well- 
understood, and many of these schemes perform competitively when integrating nonlinear evolution 
equations danang. Despite these advances, no practical high-order exponential integrators have 
been developed. High-order ETD Adams methods are largely unusable due to their small stability 
regions, and there are no ETD Runge-Kutta schemes of order greater than five. 

Nevertheless, high-order exponential integrators could prove useful if paired with spatial spectral 
discretizations, especially on periodic domains. Spectral methods exhibit exceptional accuracy and have 
been shown to be remarkably successful when applied to nonlinear wave equations mmolll]. When 
applying spectral methods on PDEs with smooth solutions, the time integrator often limits the overall 
order of accuracy. The development of stable, high-order integrators will allow for more accurate 
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numerical simulations at reduced computational costs and will better balance spatial and temporal 
accuracy. 

In order to develop high-order ETD schemes, we turn our attention to spectral deferred correction 
methods (SDC), originally developed by Dutt, Greengard, and Rokhlin [8]. SDC methods are a class of 
high-order, self-starting time integrators for solving ordinary differential equations. By pairing Euler’s 
method with a Picard integral equation, SDC methods achieve an arbitrary order of accuracy and 
favorable stability properties. Remarkably, they are simple to implement, even at high order. In the 
past decade, there has been a continuing effort to analyze and improve these methods [38l [39l [6l 

[m Hi nil Us]. In particular. Minion introduced implicit-explict spectral deferred correction schemes 
(IMEXSDC) for integrating stiff semilinear systems [31]. To date, these methods remain the only 
practical arbitrary-order IMEX integrators. 

In this paper, we present a new exponential integrator based on spectral deferred correction methods. 
Our new integrator, which we call ETD SDC, allows for an arbitrary-order of accuracy, has favorable sta¬ 
bility properties, and outperforms state-of-the-art ETD schemes when low error tolerances are required. 
In Section we provide a brief introduction to spectral deferred correction methods before deriving our 
ETDSDC method and discussing IMEXSDC. In Section we analyze and compare the stability and 
accuracy regions of these two methods. In Sectionwe discuss two techniques for accurately initializing 
the coefficients for our ETDSDC method. Einally, in Sectionwe perform numerical experiments com¬ 
paring our ETDSDC method against IMEXSDC and ETDRK4, a well-known fourth-order exponential 
integrator [7j. 


2 Spectral Deferred Correction Methods 

In this section, we provide a review of Euler-based spectral deferred correction methods |8|, before 

and the IMEXSDC method [ST] in Section |2.4| To 
introduce SDC methods, we consider a first-order initial value problem of the form 

cP\t) = F{t,cP) 

4>{a) = 4>a 

where (p E and F{t, (p) is v times differentiable for v ^ 1. We then shift our attention to a semi-linear 
first-order initial value problem of the form 

(p'it) = A(p + Af{t,(p) 

(p{a) = (pa 

where again cp e Af G C', and A is a d x d matrix (not necessarily diagonal). The continuity 
conditions on Af{t, (p) and F{t, (p) are stronger than the Lipschitz continuity required for existence and 
uniqueness, but they ensure that high-order methods can be applied successfully. 

2.1 Preliminaries 

Spectral deferred correction schemes iteratively improve the accuracy of an approximate solution to 
Eq. 0 by repeatedly solving an integral equation that governs error. This integral equation is of the 
form ^ 

y{t)=y{a)+f g{s,y{s))ds + r{t), (3) 

J a 

where r(a) = 0. As first proposed by Dutt et al. [8], we can approximate the solution to Eq. <© at 
points to, ti, ..., tm using the implicit (-t = 1) or explicit {£ = 0) Euler-like method 

y{tn+l) = y{tn) + hng{tn+e,y{tn+e)) + r{tn+l), (4) 

where hn = tn+i - tn- 


deriving our ETDSDC method in Section 2.3 
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To arrive at the error equation of the form ([^, we let be an approximate solution to Eq. 0 . 

and let the error be E{t) = (/)(t) — By considering the integral form of Eq. 0 , one arrives at 

4>{t) = 4>{a) + [ F{s,(j){s))ds. 

J a 

Substituting (p{t) = + E(t) leads to the integral equation 

E{t) = + 0'=(a) + E{a) + f E{s, + E{s))ds. (5) 

J a 


Introducing the residual 


R{t,a,(p'^) = 


(p'‘{a)+ [ F{s,4>^{s))ds 

J a 


-^\t) 


allows us to rewrite Eq. 0 


E{t)=E{a)+ f G{s,E{s))ds + R{t,a,4>'=), 
J a 


(6) 

(7) 


where 

G{s, E{s)) = F{s, <i>\s) + E{s)) - F{s, (8) 

Rewriting Eq. 0 in this manner isolates the residual and the error terms and leads to an equation of 
the form The residual R{t,a^(l)^) depends only on known quantities and can be approximated to 
arbitrary accuracy via numerical quadrature of the function F(t, (j)^{t)). If we consider a single timestep 
of method Q applied to Eq. 0. and suppose that is a sufficiently good approximation so that 


sup ||E’(t)|| = 0{h^) for h = tn+i — tn and m G N, 

tG [^Ti+l 


then, since F{t^(j)) is Lipchitz continuous in we have that 

\\hG{s, E(s))|| = h||F(s, + E{s)) - F{s, <^'=(s))|| = 0(h™+i). 

Thus, the Euler-like method Q is sufficient for estimating E{t) to in the interval [tnRn-\-i]- 

This approximate error, which we denote by E^{t), can be used to obtain an accurate solution 

— E^{t). This process can be repeated M times to obtain a sequence of increasingly 
accurate approximations to Eq. 0 . 

To implement this strategy numerically, Dutt et al. proposed to divide each timestep [tnRn-\-i] into 
N substeps or quadrature nodes which we denote via ..., [H]- This enables us to represent the 

approximate solution as an interpolating polynomial which passes through the quadrature points. 
We can then calculate a provisional solution (t) at each node using either forward or backward Euler, 
and obtain a sequence of higher-order approximations E E^~^{tnj) by repeatedly 

approximating the error E{t) at each quadrature node using 0. 

The choice of the nodes ... Rn,N affects the quality of the quadrature approximation used to 
determine Eq. 0 . Dutt et al. use Gauss-Legendre points, and Minion has studied the implications of 
using different quadrature nodes [27]. After M correction sweeps, the order of accuracy at each node is 
min(A/', M + 1), regardless of the choice of quadrature nodes [T4l l39] . 

To simplify our discussion, we consider only a single timestep of spectral deferred correction from 
tn = 0 to tn+i- We find it most convenient to describe SDC methods in terms of normalized quadrature 
points which reduce to the quadrature points if the stepsize h = 1. Throughout the rest of this paper 
we will make extensive use of the following definitions: 


Stepsize: h = t^+i — tn Normalized nodes: = tn.ij^ 

Substeps: hi = tn,i+i — tn^i Normalized substeps: r]i = hijh 


We will use the notation to denote a spectral deferred correction method which uses the 

quadrature points and performs M correction sweeps. Eor brevity we also use the variables 

= (t)^{tn^i) and Ef = E^{tn,i) to denote the approximate solution and the error at the ith quadrature 
node after k correction sweeps. 
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2.2 Euler-Based Spectral Deferred Correction Methods 

We now describe Euler-based spectral deferred correction methods in detail. Implicit and Explicit SDC 
methods use Implicit or Explicit Euler respectively to determine the provisional solution at the 

quadrature points Applying the Euler-like method @ to Eq. 0 one obtains an approximation 
of the error E{t) at each of the quadrature points. Every step of this Euler-like method requires 
approximating the residual term; we describe this process below. 

Approximating the Residual Term: During the correction sweep, is known at the quadra¬ 
ture points. The residual term 0 can be approximated for t = and a = hri at the cost of N 

function evaluations F{hTi^(j)^) via 


R{hTi+i,hTu^'') = <f>\hTi) - <t>\hTi+^) + 
where denotes the Nth order numerical quadrature approximation to 

phTi+i 

/ F{s,4>'^{s))ds. (9) 

J hri 

The coefficients for this numerical quadrature can be obtained for general quadrature points using 
an algorithm which we propose in Section® Eor Chebyshev quadrature points, a fast 0{N \og{N)) 
matrix-free algorithm exists for computing ^ [33] . 

Given the initial condition (/){ = (/)(a), we can express a single timestep of an SDCjlf method algorith¬ 
mically: 


Implicit {£ = 1) or Explicit {£ = 0) SDC]^ 

Note: E\ = 0 

• Initial Solution (Euler): 


for i=l to N-1 

= <!>} + hiF{hTi+i, 


• Correction & Update: 
for k=l to M 
for i=l to N-1 

Ei+i = Ef + hiG{hTi+e, E’^, 4 A) + E(hTj+i, hr,, cf)'") 

</>f+/=0r+i+^f+i 



By substituting the expression for into the update formula for , noting that -h Ef, 

and using Eq. 0, one arrives at the following direct update formula: 

[Fihn+e, ) - F{hTi+e, cf>'y+,)] + 7 ^'(</>''). 

This compact form for spectral deferred correction methods was first mentioned in [31] but was not 
recommended due to potential numerical rounding errors. However, in our numerical experiments, we 
find that this compact formula leads to simpler codes and equally accurate results. We therefore make 
use of this compact update formula in all of our codes. 

2.3 ETD Spectral Deferred Correction Methods 

We now introduce a new class of exponential integrators based on spectral deferred correction for solving 
Eq. 0> which we repeat here for convenience: 

= A(j) + M{t, (j)), 

4>{a) = (j>a- 
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To derive ETD spectral deferred correction schemes, we seek an error equation of the form 


y{t) = y{a)e^^* 9 {s,y{s))ds + r{t). 

J a 


(10) 


We propose to approximate the solution to Eq. ( [1Q| ) by replacing g{s^y{s)) with a one-point approxi¬ 
mation, leading to the explicit {£ = 0) or implicit (^ = 1) ETD Euler-like method 


y{tn+i) = y{tn)e^^ + A ^ - I] g{tn+e, y{tn+e)) + ritn+i). 


( 11 ) 


To arrive at an error equation of the form (10), we let be an approximate solution of Eq. (|^, and 

define the error to be E{t) = Applying variation of constants, we obtain the integral form 

ofEq. (H, 

/•* 

</.(t) = + / e^(*-*W(s,<(.(s))ds. 

J a 

Substituting 4’{t) = + E{t) leads to the integral equation 

E{t) = + E{a)) + f (p'^is) + E{s))ds. (12) 

J a 

Introducing the residual 


{s))ds 

J a 


cp\t) 


(13) 


allows us to rewrite Eq. (12) as 


where 


E{t) = E(a)e'^(*““) + f E{s))ds + Re{t, a, y), 

J a 


H{s, E{s)) = Af{s, y (s) + E{s)) - V(s, (p’‘{s)). 


(14) 


(15) 


Now that we have obtained an error equation of the form ( |1Q[ ), we are free to proceed in the same 
manner as Euler-based spectral deferred correction. The provisional solution is calculated at the 
quadrature points using either implicit or explicit ETD Euler and the error at each quadrature point is 
estimated using (11). As before, we describe the computation of the residual term. 


Approximating the Residual Term: During the correction sweep, is known at the quadra¬ 
ture points. The residual ( p^ can be approximated for t = and a = hri at the cost of N function 

evaluations via 

ReihTi^i, hTi, = (j)^{hTi)e^^^ - 0^(hri+i) + (16) 

where denotes the weighted N point numerical quadrature approximation to 





(17) 


where the weight function is w{s) = eA(^*+i We describe in detail how to obtain the coefficients for 
this weighted quadrature in Section 

We use ETDSDC]!/ to denote an ETD spectral deferred which performs M correction sweeps on the 
quadrature points {ri}fLi. Given the initial condition = 0(a), we can express a single timestep of 
an ETDSDC]!^ method algorithmically: 
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Implicit {£ = 1) or Explicit {£ = 0) ETDSDC]^ 

Note; E^ =0 

• Initial Solution (ETD Euler): 


for i=l to N-1 


(pl+i = - I] M{hTi+i,, (/))+^) 


• Correction & Update: 


for k=l to M 


for i=l to N-1 


Ef+i = + A-i - /] H{hTi+t, E\ </.'=) + R^{hTi+uhTu 



By substituting the expression for into the update formula for noting that + £’' 


A^+l - 




and using Eq. (15), one arrives at the following direct update formula: 




+ A 


-1 


hiK 


- 1] [mn+e, <p1+l) - mn+e, + Wi+\<pX- (18) 


Though we have derived both an implicit and explicit exponential integrator, we will be solely consid¬ 
ering the explicit exponential integrator throughout the rest of this paper. 


2.4 IMEX Spectral Deferred Correction 

We now briefly discuss Minion’s IMEXSDC]!^ method for solving Eq. @ m- The provisional solution 
is calculated using IMEX Euler. The error and residual equations can be derived by repeating 
the procedure outlined in Section 2.1 with F(t^y) = Ky M{t^y). This leads to 


E{t)=E{a)+ f [AE{s) + G{s,E{s))]ds + R{t,a,4>'=), 

J a 

H{s, E{s)) = V(s, E{s) + 4>'^{s)) - V(s, 4>''{s)) 
R{t,a,(l)'^)= cp'^{a)+ f [Acp'^{s)+Af{s,cp'^{s))]ds - 

J a 


Notice that Eq. (19) is of the form 

y{t)=y{a)+ f [Ay{s)+g{s,y{s))]ds + r{t). 
J a 


(19) 

( 20 ) 
( 21 ) 

( 22 ) 


We can approximate Eq. (22) by treating the linear term implicitly and the nonlinear term explicitly, 
yielding the IMEX Euler-like scheme 


y{tn+i) = {I - hA) ^ [y{tn) + hg{tn+i, y{tn+i)) + r{tn+i)] . 

The residual term ( [21] ) is approximated exactly as described in Section |2.2[ except the integrand in 
Eq. ^ is now Ac/)^(s)N{s^ (s)) . We denote the quadrature approximation to the residual for 
IMEXSDC by R{t^a^(j)). Given the initial condition (j)\ = 0(a), we can express a single timestep of an 
IMEXSDC]!^ method algorithmically: 
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IMEXSDCj)? Method 

Note; Ef = 0 

• Initial Solution (IMEX Euler): 

for i=l to N-1 

<^1+1 = [/ - hiA]-^ + hiNihn, <pl)] 


• Correction & Update: 
for k=l to M 
for i=l to N-1 

El, = [I- hiA]-^ \E'y + hiH{hTu E, + R{hTi+,, hn, 4>l 

4>Ut=^i+i+Eli 



By rewriting the error formula implicitly so that 

Ef+i = [e^ + + H{hTu E, </.'')) + R{hTi+^,hTi, 4>’‘) 

substituting this expression into the update formula for (j)i+iJ S'lid noting that 

El, = - cf>l„ cf>^+^ = + E^ 

one arrives at the following direct update formula: 

= [/ - hiA]-^ [./.yi - {hiA)cf>i, + hi{M{hTi, - M{hTi, 0f)) + ii+\cf>i 
where denotes the numerical quadrature approximation to 


rhTi- 
J hri 


K(j)^{s) ^ J\f{s,(j)^{s))ds. 


3 Stability and Accuracy 

Determining the stability properties of IMEX and ETD integrators is non-trivial. A commonly used 
approach is to consider the model problem 


0' = /i0 + \(j) 

m = 1 


(23) 


where /i, A G C and the terms X(j) act as the linear and nonlinear term respectively. This model 
problem highlights stability for Eq. <§ when it is possible to simultaneously diagonalize both the linear 
and nonlinear operators around a fixed point. Though this analysis does not extend to general linear 
systems, it has proven useful for predicting stability properties of IMEX and ETD methods on a variety 
of partial differential equations 


Applying an ETDSDCj!^ or IMEXSDC]!/ method on Eq. (23) leads to a recursion relation of the 
form 

<t){tn+l) = IpN (r, z)<Pitn) 


where r = fih^ z = \h^ and h denotes the timestep. As with all one-step methods, the stability region 
is defined as 

S = {{r,z) e C^ \ipl{r,z)\ < 1}. 

We list the stability functions (r, z) for ETDSDCjif and IMEXSDC])f schemes in Table 

We choose to analyze stability for PDEs with linear dispersion and dissipation; thus, r = h/i and 
z = hX are complex-valued. Several strategies have been proposed for effectively visualizing the resulting 
four-dimensional stability region. As in 0171 US], we choose to overlay two-dimensional slices of the 
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ETDSDC Stability Functions 


^Hr,z) = 1 


prr]i _ 1 

'ipj+i = +- - - 

N 

prr]i _ 1 

i’t+i = ^ 

^ i=i 

where = f e^^'^^+^-^^Lj{s)ds, Lj{s) = JJ 

Jri i^i \^j ~ ^l) 

W 


(24) 


IMEXSDC Stability Functions 


tl)^{r,z) = 1 




1 

i+l 


/c+1 

i+1 


where 


^ 1 + 

I-rr]i 




'^i + “ '^i) ~ f'di'^i +1 + {'r + z) 'iiji’j 

1 -rrji 



N 

L,(s)=n 


1=1 


i^j 


(■s - n) 

-n)' 


(25) 


Table 1: Stability functions for ETDSDCjy and IMEXSDCjy methods. As r —J- 0 the stability functions 
of both methods limit to that of an explicit SDCjif method. 
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stability regions, each corresponding to a fixed r value. For the sake of brevity, we focus our attention 
on 8th order methods where N = 8, M = 7 and on 16th order methods where N = 16^ M = 15. For 
all methods, we select the Chebyshev quadrature nodes 


1 q f -1) 
’■• = 2(,*““n"Frr 


i = 1 , ..., Af. 


We pick a range of real, imaginary, and complex r values to simulate nonlinear PDFs with varying 
degrees of linear dispersion and dissipation. We plot stability regions pertaining to 


r G —1 • [0,30], r G li • [0,30], and r G exp(37ri/4) • [0,30] (26) 

in Figure For these three r ranges, we find that the stability regions of all methods grow as |r| 
increases. For imaginary r, the stability regions for ETDSDC methods temporarily decrease before 
growing. Though all methods exhibit satisfactory stability properties, IMEXSDC methods allow for 
coarser timesteps on a wider range of (r, z). Overall, our results suggest that both IMEXSDC methods 
and ETDSDC methods exhibit good stability properties on a wide range of stiff nonlinear evolution 
equations. 

When analyzing spectral deferred correction methods, it is also common to plot accuracy regions. 
Accuracy regions highlight the restrictions on the stepsize h so that error after one timestep is smaller 
than e > 0. They are simply defined as 

A = {{r,z) e C^, \tppf{r,z) - exp(r + z)\ < e}. 


They were introduced in [8] for comparing the efficiency of high-order methods, and provide a more 
detailed picture than stability regions which solely differentiate between convergent and divergent (r, z) 
pairs. 

We find that as |r| increases, the accuracy region containing z = 0 decreases rapidly for ETDSDC]!^ 
methods and vanishes entirely for IMEXSDC]!^ methods. This behavior can be understood from 
Eq. (24) and Eq. (25). Eor the ETDSDCjlf methods it follows that (r, 0) = exp(rh); moreover, 
since the stability function (r, z) is continuous, then for any e > 0, there exists a nontrivial accuracy 
region surrounding z = 0. The same cannot be said for IMEXSDC schemes since Eq. ( [^ satisfies the 
weaker relation '0(r, 0) = exp(rh) + 0{rh); hence, as r becomes sufficiently large, there need not exist 
an accuracy region around z = 0. 

We present accuracy regions for e = I x 10“^ in Figure We consider the three ranges of r values 
in ( [ 2 ^ , but due to rapidly shrinking accuracy regions, we are only able to visualize different subsets 
of r values for each numerical method. ETDSDC]!^ schemes outperform IMEXSDC]!^ schemes for all 
tested values. Accuracy regions for the ETD methods decrease more slowly, and the non-vanishing 
accuracy regions around z = 0 guarantee accuracy for any r so long as z is chosen sufficiently small. 
The MATLAB code used to generate these figures can be found in [5] and can be easily modified to 
generate stability and accuracy plots for other ETDSDC or IMEXSDC methods. 


4 Calculating 

Every iteration of an ETDSDC]!^ method requires computing which denotes the weighted 

quadrature approximation to 

nhrij^i 

J hTi 

To arrive at a formula for we let N/((/)) = M{hTi^(j){hTi)) and replace N{s^(j)^{s)) in 

Eq. ^ with the Lagrange interpolating polynomial L{s) that passes through the quadrature points 
so that 


rhTi^i ^ 

JhTi 


(27) 
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Stability Region Plots 


’ = J^o/2 


r = 2R 


0 


ETDSDC| 


IMEXSDC| 


ETDSDC}! 


IMEXSDC}! 



Figure 1: Stability regions for Sth order and 16th order methods with Chebyshev quadrature nodes. 
Colored contours correspond to different r values as described in the legend. We plot an additional 
black contour for the ETDSDC jg method on the dispersive model problem to show that stability regions 
eventually grow for sufficiently large imaginary r. For large |r|, increasing the order of the ETD and 
IMEX methods does not lead to significantly larger stability regions. 
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Accuracy Region Plots 




ETDSDC| 


IMEXSDC| 


ETDSDC}! 


IMEXSDC];| 


Dissipative Model Problem: r € — 1 • [0,30] 



Dispersive Model Problem: r € li ■ [0, 30] 


Ho = 

(M 

II 

0 

Rq = 15i 

2 

"16 

16 




■ 


8 

0 


0 




.1 


-8 

0 



, -2 

._ ^^^^ -16 

, -16 


-2 -1 0 1 2 -2 -1 0 1 2 -16 -8 0 8 16 -16 


Ho = 9z 



-8 0 8 


16 


Dissipative/Dispersive Model Problem: r G exp(37ri/4) • [0,30] 


Rq = 


Rq = 





Figure 2: Accuracy regions corresponding to e = 1 x 10“^ for 8th order and 16th order methods with 
Chebyshev quadrature nodes. Colored contours correspond to different r values as described in legend. 
We choose Rq in each figure so that the red contour marks a near vanishing accuracy region around 
2 ; = 0. As expected, 16th order methods possess larger accuracy regions for a wider range of r than 8th 
order methods. 
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For low-order methods, explicit formulae for Wi^i can be derived by forming L{s) and repeatedly applying 
integration by parts. Unfortunately, this direct calculation leads to increasingly involved formulae for 
large N. We therefore seek a general procedure for determining Wi^i for any N. We propose to express 
the weights Wi^i in terms of the well-known functions 


<Pn{z) = 7—^TT [ ^d(j. 

- 1 )! Jo 


in- 

using a stable algorithm developed by Fornberg for determining finite difference coefficients [9]. We 
describe our algorithm in Section |4.1[ before discussing (f functions and two well-known methods for 


initializing them in Section 4.2 


4.1 Proposed Algorithm 

To arrive at a convenient expression for we propose to apply the change of variables 


s = h [(ri+i - Ti)(j + Ti] = hi(T + hTi, 


(28) 


to the integral term in (27), expand the Lagrange interpolating polynomial L{s{a)) as a Taylor poly¬ 
nomial, and rewrite the result in terms of ip functions. Applying the change of variables (28) leads 
to 1 1 

hi f L{s{a))da = hi f Pi{cF)dcF, 

Jo Jo 

where Pi{(j) is the Lagrange interpolating polynomial which passes through the points 


{{qi,h^i{(t)^))}^=i and Qi^i = {n -ri)/{ri^i - Ti) 

denote the scaled, translated quadrature nodes hri under the transformation ( [^ . Next, we define the 
finite difference coefficients so that 




N 


Expanding Pi{cr) as a Taylor polynomial we obtain 

j=0 1=1 


N-i r N 


da. 


Reordering terms we arrive at 

= hi 


N 


1=1 

N 


N-1 




= hi 


1 = 1 


j=0 

N-1 




hL f phiAil- 

n Jo 


^^a^ da 


[aJp^i+^ihiK) 


j=0 


By defining the functions 

N-1 

Wi,iiz) = hi ^ afj(pj+i{z), 

j=0 

we obtain a convenient expression for the weighted quadrature rule: 

N 




1 = 1 


( 29 ) 
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To successfully implement this procedure, we must determine the finite difference coefficients a/- [ and 

the matrix functions (pn{hiA). The coefficients a^- \ can be rapidly obtained using the stable algorithm 
presented in [9]. We define the functions: 

• weights{zo^ [qi,... ^qn],m): returns a finite difference matrix a for computing m derivatives at zq, 
assuming qj are the quadrature points. This calling sequence is consistent with the implementation 

in [To]. 

• initPhi{z^n)\ returns the functions (pi{z) for i = 0,...,n. We discuss two possible implementa¬ 
tions in Section [4^ 

The algorithm for computing Wi^i{z) for an ETDSDC]!/ method can be written as: 


Computing Wi^i{hiA) 
for i=l to N 

[ipo{hiA), ..., ipN{hiA)] = initPhi(hiA, N) 
for j=l to N 

= {Tj -Ti)l{Ti+l - n) 
a(*) = weights(0, [qi,... qn], N - 1) 
for 1=1 to N 
for j=0 to N-1 

Wi^i{hiA) = Wi^i{hiA) + a^jipjj^i{hiA) 


When computing Wi^i{hiA), it is convenient to save ipo{hiA) and ipi{hiA) since both are required for 
the ETD Euler method. 


4.2 (p Functions 

The coefficients of all exponential integrators can be expressed in terms of p functions [snuaisniEi]- 
The nth (p function can be defined in the following ways; 


(30) 

(31) 

(32) 

The first few Pn{^) are given by 

_ 1 p^ _ 1 _ p^ — 1 — 2^ — — 7^ 

Mz) = e^, Mz) = -, Mz) = -^—, miz) = -^ 

We can now rewrite the compact update formula ( [T^ as 

+ MhiA) [N{hn+e,<l)pl) - iV(hn+,,((.f+,)] + 


Integral Form: 


Series Form: 


Reeursion Relation: 


‘Pniz) = < 


n = 0 

' , ^ n>0 

(n - 1)! Jo 

OO /j. 

" £ {kT^. 


k=0 


( \ (n-1)! . X 2 


Erom their series definition, it follows that the functions Pn{z) are entire; nevertheless, it is well-known 
that explicit formula for Pn{^) are prone to catastrophic numerical roundoff error for small jzj. Various 
strategies for overcoming this difficulty have been compared extensively [T] . We briefly outline a method 
based on scaling and squaring m and a method based on contour integration [22]. Other approaches 
involve Krylov subspace approximations mi US] and improved contour integrals im but we do not 
consider them in this paper. 
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4.2.1 Taylor/Pade Scaling and Squaring Algorithm 


The scaling and squaring algorithm for calculating ip functions is a generalization of a well-known 
algorithm for computing matrix exponentials US]. For small |z|, Pn{z) can be accurately evaluated via 
the Taylor series (31) or via the diagonal (m,m) Fade approximation, whose explicit formula is given 
in [37]. This initial approximation can be used to obtain Pn{z) for large | 2 :| by repeatedly applying the 
well-known scaling relation 


^n{z) 


1 

2n 


<^0 (i) (f) + y] 

i=l 


(n — i)\ 


(33) 


We present pseudocode for an m-term Taylor series procedure for initializing <pi (A) in Table A 
MATLAB implementation of the Fade scaling and squaring algorithm is freely available in [2] and can 
be easily used to initialize Pn{^) foi* both scalar and matrix A. 


4.2.2 Contour Integration Algorithm 

An alternative algorithm for initializing ETD coefficients was first suggested in [22] . Since the functions 
Pn{z) are entire, Cauchy’s integral formula can be used to obtain p{z) at problematic regions near 
2 : = 0. We highlight this procedure for both scalar and matrix A in Table [^ assuming that the explicit 
formula for Pn{z) is known. If this is not the case, then it is convenient to combine Eq. (32) with the 
discretized contour integral so that 


(^) 


1 C' + J’e*®) - l/(^^ - 1)! 

P ^ A + 

j=0 


(34) 


This allows one to progressively evaluate Pn{^) for n = 1,..., A. Eor scalar A we use Eq. (34) when 
|A| <1 and Eq. (32) when |A| > 1. Eor matrix A we find that the technique based on scaling and 
squaring is faster and more accurate, especially for matrices with large norm. 


5 Numerical Experiments 

In this section, we numerically solve four partial differential equations in order to compare ETDSDC]!^ 
and IMEXSDC]!f methods of orders 4, 8,16, 32 against the fourth-order Runge-Kutta method (ET- 
DRK4) developed in [7]. We have chosen to include ETDRK4 in our tests since it was shown to 
perform competitively nang, and provides a good reference for comparing SDC based schemes to ex¬ 
isting ETD and IMEX methods. We provide our MATLAB and Eortran implementation of ETDSDCjlf, 
IMEXSDC]!^ and ETDRK4 in [5] along with code for reproducing our numerical experiments. 

In all our numerical experiments, we apply a fine spectral spatial discretization so that error is 
primarily due to the time integrator. In our first three experiments we impose periodic boundary 
conditions and solve the FDEs in Eourier space. This is convenient since it leads to an evolution 
equation of the form ([^ where the matrix A is diagonal. In our final experiment we consider a more 
challenging example where A is a dense matrix. We base our first three numerical experiments from 
[mm] so that our results can be compared with those obtained using other IMEX and ETD schemes. 

Since we consider methods of varying order, our experiments are based on the number of function 
evaluations rather than the step size h. We compute reference solutions by using four times as many 
function evaluations as used in the experiment. To avoid biased results, we average the solutions of at 
least two convergent methods when forming our reference solutions. Eor each FDE, we present plots of 
relative error vs. function evaluations, relative error vs. stepsize, and relative error vs. computational 
time, where the relative error between two solution vectors x and y is ||x — y||oo/||x||oo- Though we 
solve equations in Eourier space, we compute relative errors in physical space. We do not count the time 
required to initialize ETD coefficients in our time plots. We also make no specific efforts to optimize 
our code, thus timing results only serve as an indication and may vary under different implementations. 
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m-Term Taylor Scaling Squaring Procedure for Matrix/Scalar A 

• Select Scaling Factor: 

Let s G N so that A/2^ oo < S{m) See Appendix [a| or [24] for choosing 5{m). 

• Initialize ipi{K/2^) via Horner^s Method: 

for i=0 to N 

p. — A 1 I 

* (m+i)! ' (m+i — 1)! 

for k=0 to m-2 

Pi - APi + 

^i(A/2*) = Pi 

• Obtain (pi{A) via Eq. 

for i=l to s 
for n=0 to N 

^n{A/r-n = ^ <Po (A/2-^+1) (A/2-^+1) + x: 

Contour Integral for Scalar A < 1 

Contour Integral for Matrix A 

• Cauchy Integral Formula: 

^„(A)= ^ 

^ 27ri/r(z-A) 

• Choosing T: 

Let r = Re^^ + A for G [0, 27r]. The radius 

R should be chosen so that contour never 
comes near the origin. 

1 p‘2tt 

‘/^n(A) = 7 r <^n(-Re*® + A)de 

Jo 

• Discretization via Trapezoidal Rule: 

Let Oj = 2112 jP^ fhen for P sufficiently large, 
(Pn{A) is approximately 

p X] + RF^^). 

j=0 

For scalar A > 1, use Eq. (32). 

• Cauchy Integral Formula: 

¥’n(A) = ^ / Pn{z)izl - A)~^dz 

• Choosing T: 

Let r = Re^^ + zq for 0 G [0,27r]. The ra¬ 
dius R and center zq must be chosen so that 
contour encloses spectrum of A. 

1 

¥’n(A) = 7r Pn{Re''^ + A)d9 

Jo 

• Discretization via Trapezoidal Rule: 

Let Oj = 271 jIP, 7 j = Reyip{0j) + zq, then for 
P sufficiently large, ^n{A) is approximately 

j=0 ^ / 

where (Pnijj) initialized like scalar A. 


Table 2: Scaling & squaring, and contour integral methodology for initializing (/^^(A). 
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The results presented in this paper have been run on a 3.5 Ghz Intel i7 Processor using our double 
precision Fortran implementation. We describe each of the four problems below. 

The Kuramoto-Sivashinsky (KS) equation models reaction-diffusion systems [26]. As originally 
presented in [22], we consider the KS equation with periodic boundary conditions: 

~ '^xxxx 2 ^x ^ 

u{x, t = 0) = COS (^) (l + sin (^)) , x G [0, 647r]. 

We numerically integrate Eq. ( [3^ using a 1024 point Fourier spectral discretization in x and run the 
simulation out to t = 60. The KS equation has a dispersive linear term A with eigenvalues given by 
X{k) = where k denotes the Fourier wavenumber. We present our numerical results in Figure 

[31 


The Nikolaevskiy equation was originally developed for studying seismic waves [32] and now serves as 
a model for pattern formation in a variety of systems [36]. As originally presented in m, we consider 
the Nikolaevskiy equation with periodic boundary conditions: 

Ut = adlu + Pdlu -dl{r-{l+ dlf) u-\ {u^)^ , (36) 

u{x^ t = 0) = sin(x) + esin(x/25), x G [—TStt, TStt] 

where r = l/4, o = 2.1,;d = 0.77, and e = 1/10. We solve the Nikolaevskiy equation using a 4096 point 
Fourier spectral discretization in x and run the simulation out to t = 50. The Nikolaevskiy equation has 
a dissipative and dispersive linear term with eigenvalues given by X{k) = k‘^{r — {1 — k^^Y) — iak^ + i/Sk^^ 
where k denotes the Fourier wavenumber. We present our numerical results in Figure [^ 

The quasigeostrophic (QG) equations model a variety of atmospheric and oceanic phenomena [35] . 
As originally presented in [12], we consider the barotropic QG equation on a /d-plane with linear Ekman 
drag and hyperviscous diffusion of momentum with periodic boundary conditions, 

= - [pdxi^ + eVV + + u • V(VV)] (37) 

ip{x,y,t = 0) = ^exp (^-8 {2y^ + x^j2 - 7r/4)^) , 

{x,y) e [-7r,7r] 


where 2 p{x,y) is the stream function for two-dimensional velocity u = {—dy^p, dx^p)^ e = 1/100, and 
u = 10“^^. We run the simulation to time t = 5 using a 256 x 256 point Fourier discritization. We 
consider a different initial condition than the one presented in [12], since \/‘^'ip{x,y) was originally 
chosen to be discontinuous at the point (0,0). We note that Eq. (37) describes the change in the 
vorticity uo = in terms of the stream function pj. In order to obtain pj at each timestep, it is 

necessary to solve Poisson’s equation = uj. Since we are solving in Fourier space, it follows that 


^k,i 



k = l = 0 
otherwise 


where k and I are the Fourier wave numbers and t/, uj denote the discrete Fourier transforms of ip and 
UJ. The QG equation has a linear term with strong dissipation and mild dispersion with eigenvalues 
given by X{k^ 1) = - iy{k^ + /^). We present our numerical results in Figure 4 
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Performance Results for Kuramoto-Sivashinsky Equation 



Accuracy vs Function Evaluations 



Accuracy vs Stepsize 



10 ® 10 ® 10 ^ 10 ° 10 ^ 
step Size h 



Performance Results for Nikolaevskiy Equation 




Figure 3: Performance results for the Kuramoto-Sivianshi and Nikolaevskiy equations. Gray dashed 
lines of increasing steepness in the accuracy vs stepsize plots correspond to O(h^), 0 {h^) and 
respectively. IMEXSDC schemes experience significant order reduction on both problems. 


17 



























Performance Results for Quasigeostrophic Equation 



Accuracy vs Function Evaluations 



Accuracy vs Stepsize 

10 “ - 


10 ““ 

m 



Step Size h 


Accuracy vs Computer Time 


10 ° 

10 ““ 



Computer Time (sec) 


Performance Results for Korteweg-de Vries Equation 



Accuracy vs Stepsize 



Accuracy vs Computer Time 



Figure 4: Performance results for the Quasigeostrophic and Korteweg-de Vries equations. Dashed lines 
of increasing steepness in the accuracy vs stepsize plots correspond to O(h^), 0 {h^) 0{h^^) and 0 {h^^), 
respectively. Notice that high-order IMEXSDC schemes are unstable on the KDV equation. Order 
reduction does not occur for any method on the quasigeostrophic equation, but affects both IMEXSDC 
and ETDSDC schemes on the KDV equation. 
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The Korteweg-de Vries (KDV) equation describes weakly nonlinear shallow water waves. In 1965 
Kruskal and Zabusky observed that smooth initial conditions could give rise to soliton solutions [42]. 
As in their original numerical experiment, we consider the KDV equation on a periodic domain 

T 2^^ )^] 

u{x, t = 0) = cos(7rx), X G [0, 2] 


where 6 = 0.022 and the simulation is run out to time t = 3.6 /tt. The eigenvalues of the linear terms are 
given by X{k,l) = 6ik^; thus this equation possess a purely dispersive linear term. Unlike our previous 
examples, we solve this PDE in physical space where the resulting differentiation matrix is no longer 
diagonal. The nondiagonal case is more challenging since the coefficients Wi^n in Eq- (27) are now matrix 
functions. In practice it would be more efficient to consider a lower-order spatial description and apply 
Krylov space or contour integral techniques that avoid explicitly initializing the requisite ETD matrices. 
Nevertheless, we consider this example to test the robustness of the scaling and squaring algorithm. 
Eor IMEXSDCjlf schemes it is necessary to repeatedly solve the system Kx = / at each timestep. We 
perform an initial LU factorization of A to expedite this process. We present our numerical results for 
the KDV equation in Eigurej^ 


5.1 Discussion 

Our results demonstrate that high-order methods can lead to significant speedup when solving nonlinear 
wave equations to high accuracy. Methods of order 8 and 16 were able to achieve the smallest error 
using the fewest function evaluations and the least overall CPU time. Interestingly, the error threshold 
separating good and bad performance for high and low order methods varied significantly in each 
experiment. Overall, ETDSDC methods consistently achieved better accuracy than corresponding 
IMEXSDC methods, and did not suffer from crippling order reduction on any of the problems we tested. 
Amongst the fourth order methods, ETDRK4 is more efficient than either ETDSDC 4 or IMEXSDC 4 . 
Moreover, ETDRK4 is the fastest method for computing solutions if error tolerances are large. Methods 
of order 32 were generally less competitive than those of 8 th or 16th order, and should only be considered 
in situations where extreme precision is necessary and quad/arbitrary-precision arithmetic allow for 
relative errors significantly below 1 x 10“^^. Einally, for diagonal A, the time required to initialize 
the ETD coefficients was insignificant as compared to overall computational time even for 32nd order 
method. 

High-order ETDSDC methods continued to perform well even in the non-diagonal case, and we found 
no evidence of catastrophic roundoff error when forming the ETD matrix coefficients Wi^i{hiA). Eor 
nondiagonal A, high-order ETDSDCjlf schemes require large amounts of memory and time to initialize 
and store the N‘^ — N requisite matrices. Moreover, the expensive matrix multiplications at each 
timestep reduced their overall competitiveness. To improve the performance of ETDSDC schemes on 
higher dimensional problems with non-diagonal linear operators, it becomes essential to use techniques 
that avoid explicitly storing the ETD matrices. 

High-order IMEXSDC schemes were unstable when solving the KDV equation on fine grids in both 
physical and Eourier space. Through additional numerical testing we find that IMEXSDC schemes can 
be unstable when integrating other nonlinear wave equations with dispersive linear terms such as the 
nonlinear Schrodinger equation. 

We make several additional comments regarding our numerical experiments. The benefits of using 
high-order methods is greatly reduced if the initial conditions are not smooth, though in certain situ¬ 
ations we found that high-order methods are rendered no less efficient than lower-order counterparts. 
The size of the integration window also affects the difference in performance between high and low-order 
methods, with the high-order methods generally benefiting on larger time domains. Chaotic equations 
can cause additional complications, as small perturbations due to rounding errors grow exponentially 
and contaminate overall accuracy. This was the case for the the KS equation where we were not able 
to integrate further without damaging the quality of the reference solution. 
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6 Conclusion 


We have demonstrated that high-order ETD spectral deferred correction schemes possess excellent ac¬ 
curacy/stability properties and outperform existing ETD and IMEX methods when solving nonlinear 
wave equations to high accuracy. Our proposed methodology for initializing ETD coefficients is robust 
and can be successfully applied to ETDSDC schemes up to 32 order accuracy, even for equations with 
non-diagonal linear operator A. We have also highlighted the advantages of ETD spectral deferred cor¬ 
rection methods as compared with IMEXSDC schemes. Our new ETD schemes consistently outperform 
their IMEX counterparts, do not appear to suffer from crippling order reduction, and retain stability 
on equations with dispersive linear terms. # 
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A Choosing 5{m) 

We describe a simple choice for S{m) from Table more sophisticated alternatives are developed in [24]. 
Let A be a matrix or scalar and let A = A/2^. We seek an integer s so that (^n(A) can be initialized 
via its mth order Taylor series without admitting an error larger than e, assuming exact arithmetic, 
can be approximated by 


^n(A) 


m 


E 

k=0 


A/e 

{k + n)\ 


+o(iiAr+i). 


The error En{A) = ||(/?„(A) — (/?™(A)|| can be expressed as 


En{A) 


E 

k=m-\-l 


A* 

{k + n)\ 


OO 

^ E 

k=m-\-l 


IIAf 

{k + n)\' 


Eor any n, En{A) can be bounded above by 

^ ||Ar+’» _ ||A|r+iexp(||A|y 

Assuming exact arithmetic, we can guarantee that En{A) < e for any s so that A = A/2^ satisfies 

^ _ ||A|r+iexp(||A||) _ 

(m + 1)! 

To avoid solving a nonlinear system for each m and e, we fix m = 20, e = 1 x 10 and let ||A|| < p. 
Solving the cooresponding nonlinear equation for p, leads to p = lA ^ 1.0. Therefore our condition on 
A reduces to ||A|| = ||A/2^|| < 1. In our numerical codes we choose ||A|| = max{||A||i, ||A||oo}- This 
leads to the condition 

ln(max{||A||i,||A||oo}) 

s - max |0, 
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